#include <iostream>
#include <vector>
#include <cmath>
// -------------------- OpenMesh
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>

struct MyTraits : public OpenMesh::DefaultTraits
{	
	template <class Base, class Refs> struct VertexT : public Base
	{
		double distance; 
	};
};
typedef OpenMesh::TriMesh_ArrayKernelT<MyTraits>  MyMesh;

bool on_edge (MyMesh::Point & p1, MyMesh::Point & p2, double value1, double value2, double isovalue, MyMesh::Point & p)
{
	std::cout<<value1 << ' '<<value2<<std::endl;
/*
	if (std::abs (isovalue - value1) < 1e-7) 
	{
		p[0] = p1[0]; p[1] = p1[1]; p[2] = p1[2];
		return true;
	}
	if (std::abs (isovalue - value2) < 1e-7) 
	{
		p[0] = p2[0]; p[1] = p2[1]; p[2] = p2[2];
		return true;
	}
*/
	if (isovalue > value1 && isovalue < value2)
	{
		p[0] = p1[0] + (p2[0] - p1[0])* (isovalue - value1)/ (value2-value1); 
		p[1] = p1[1] + (p2[1] - p1[1])* (isovalue - value1)/ (value2-value1);
		p[2] = p1[2] + (p2[2] - p1[2])* (isovalue - value1)/ (value2-value1);
//		std::cout<<"Cut1 : "<< p[0]<<' '<<p[1]<<' '<<p[2]<<std::endl;
		return true;
	}
	else if (isovalue < value1 && isovalue > value2)
	{
		
		p[0] = p2[0] + (p1[0] - p2[0]) * (isovalue - value2)/ (value1 - value2); 
		p[1] = p2[1] + (p1[1] - p2[1]) * (isovalue - value2)/ (value1 - value2);
		p[2] = p2[2] + (p1[2] - p2[2]) * (isovalue - value2)/ (value1 - value2);
//		std::cout<<std::endl;
//		std::cout<<p1[0]<< ' '<<p1[1]<<' '<<p1[2]<<std::endl;
//		std::cout<<p2[0]<< ' '<<p2[1]<<' '<<p2[2]<<std::endl;
//		std::cout<<"Cut2: "<<p[0]<<' '<<p[1]<<' '<<p[2]<<std::endl;
		return true;
	}
	return false;
}
int main(int argc, char **argv)
{
  MyMesh  mesh;


  // check command line options


  // read mesh from stdin
  if ( ! OpenMesh::IO::read_mesh(mesh, "in.off") )
  {
    std::cerr << "Error: Cannot read mesh from "  << std::endl;
    return 1;
  }


  // this vector stores the computed centers of gravity
  std::vector<MyMesh::Point>  v_result;
  std::vector<MyMesh::Point>::iterator cog_it;


  // smoothing mesh argv[1] times
  MyMesh::VertexIter          v_it, v_end(mesh.vertices_end());
  MyMesh::VertexVertexIter    vv_it;
  MyMesh::Point               cog;
  MyMesh::Scalar              valence;

  MyMesh::FaceIter   f_it, f_end(mesh.faces_end());
  MyMesh::FaceVertexIter fv_it,fv_it_end;


  int c = 0;
  /*Computing the geodesics on plane by Euclidean distance simplyly*/
    for (v_it=mesh.vertices_begin(); v_it!=v_end; ++v_it)
	{
		cog = mesh.point(v_it);
		//double d = v_it->get_distance();
		mesh.data(v_it).distance = std::sqrt (cog[0] * cog[0] + cog[1]* cog[1] + cog[2]*cog[2]);
		std::cout<<c++<<':'<< cog[0]<<' '<<cog[1]<< ' '<<cog[2]<<' '<<mesh.data(v_it).distance <<std::endl;
	}

	double isovalue = 3;
 
  for (f_it = mesh.faces_begin(); f_it!=f_end; ++f_it)
  {
	  fv_it = mesh.fv_iter(f_it);
	  fv_it_end = fv_it;
	  int count = 0;
	  MyMesh::Point triangle[3];
	  double dis [3];
	  do
	  {
		  cog = mesh.point(fv_it);
		  //std::cout<<'('<<cog[0]<<' '<<cog[1]<< ' '<<cog[2]<<')';
		  triangle[count][0] = cog[0]; triangle[count][1] = cog[1]; triangle[count][2] = cog[2];
		  dis [count] = mesh.data(fv_it).distance; 
		  if (++count >= 3) break;
		  ++fv_it;
	  }while (fv_it !=fv_it_end);
	  //std::cout<<std::endl;
	  
	  std::cout<<"Triangle :"<<std::endl;
	  for (int k = 0; k < 3; ++k)
	  {
			  std::cout<<triangle[k][0]<<' '<<triangle[k][1]<<' '<<triangle[k][2]<<' '<<dis[k]<<std::endl;
	  }
	  MyMesh::Point p;
	  count = 0;
	  MyMesh::Point linep[3];

	  if (std::abs (isovalue - dis[0]) < 1e-7)
	  {
		  linep[0] = triangle [0];
		  count ++;
		  if (on_edge (triangle[1], triangle[2], dis[1], dis[2], isovalue, p ))
		  {
			  linep[1] = p; count++;
		  }
	  }
	  else if (std::abs (isovalue - dis[1]) < 1e-7)
	  {
		  linep[0] = triangle [1];
		  count ++;
		  if (on_edge (triangle[0], triangle[2], dis[0], dis[2], isovalue, p ))
		  {
			  linep[1] = p; count++;
		  }
	  } 
	  else if (std::abs (isovalue - dis[2]) < 1e-7)
	  {
		  linep[0] = triangle [2];
		  count ++;
		  if (on_edge (triangle[1], triangle[0], dis[1], dis[0], isovalue, p ))
		  {
			  linep[1] = p; count++;
		  }
	  } 
	  else
	  {
		  if (on_edge (triangle[0], triangle[1], dis[0], dis[1], isovalue, p ))	  
		  {
			  linep[count][0] = p[0];linep[count][1] = p[1];linep[count][2] = p[2];
			  count ++;
		  }
		  if (on_edge (triangle[1], triangle[2], dis[1], dis[2], isovalue, p ))	  
		  {
			  linep[count][0] = p[0];linep[count][1] = p[1];linep[count][2] = p[2];
			  count ++;
		  }
		  if (on_edge (triangle[2], triangle[0], dis[2], dis[0], isovalue, p ))	  
		  {
			  linep[count][0] = p[0];linep[count][1] = p[1];linep[count][2] = p[2];
			  count ++;
		  }
	  }
	  if (count == 2)
	  {
		  v_result.push_back(linep[0]); 
		  v_result.push_back(linep[1]);
	  }
	  std::cout<<"count : "<<count<<std::endl<<std::endl;;
  }

  std::cout<<"OFF"<<std::endl;
  std::cout<< v_result.size()<< ' '<<v_result.size() /2 <<' '<<0<<std::endl;
  for (unsigned i = 0; i < v_result.size(); ++i)
  {
	  std::cout<<v_result[i][0]<<' '<<v_result[i][1]<<' '<<v_result[i][2]<<std::endl;
  }

  for (int i = 0 ; i < v_result.size() ; i+=2)
  {
	  std::cout<<2<<' '<<i<<' '<<i+1<<std::endl;
  }

      //if ( !mesh.is_boundary( v_it ) )
      //  mesh.set_point( v_it, *cog_it );


  // write mesh to stdout
//  if ( ! OpenMesh::IO::write_mesh(mesh, argv[3]) )
//  {
//    std::cerr << "Error: cannot write mesh to " << argv[3] << std::endl;
//    return 1;
//  }

  return 0;
}
